Random walk diffusion simulations in semi-permeable layered media with varying diffusivity

In this paper we present random walk based solutions to diffusion in semi-permeable layered media with varying diffusivity. We propose a novel transit model for solving the interaction of random walkers with a membrane. This hybrid model is based on treating the membrane permeability and the step change in diffusion coefficient as two interactions separated by an infinitesimally small layer. By conducting an extensive analytical flux analysis, the performance of our hybrid model is compared with a commonly used membrane transit model (reference model). Numerical simulations demonstrate the limitations of the reference model in dealing with step changes in diffusivity and show the capability of the hybrid model to overcome this limitation and to offer substantial gains in computational efficiency. The suitability of both random walk transit models for the application to simulations of diffusion tensor cardiovascular magnetic resonance (DT-CMR) imaging is assessed in a histology-based domain relevant to DT-CMR. In order to demonstrate the usefulness of the new hybrid model for other possible applications, we also consider a larger range of permeabilities beyond those commonly found in biological tissues.

. A grid of walker positions x near a permeable barrier located at x b . The concentration at x p , i.e. U(x p ,t + 2δt), is composed of the contribution of three different walker positions at time t (two time steps prior) through four different paths. The green path reflects at the barrier with probability 1 − p t,L→R , while the orange path passes through the barrier with transit probability p t,R→L .
Consider the probability of finding a walker at x p as the sum of probabilities of neighbouring walkers jumping to x p . Away from interfaces, this results in U(x p ,t + δt) = 1 2 U(x p − δ x,t) + 1 2 U(x p + δ x,t) .
Selecting x p to be near a barrier at x b (see figure 9 for the possible paths of walkers to reach this position) results in additional terms involving p t (both p t,L→R and p t,R→L ). The barrier location is chosen to coincide with a grid point. A walker located here will proceed left/right according to the appropriate probability of transit. The contributions of different U(x,t) terms to U(x p ,t + 2δt) in figure 9 (with the path colour indicated) are: One may now express x p relative to x b , e.g. the contribution of the orange path from the right side of the barrier is: This allows us to perform a Taylor series expansion of every term around x b , taking care to expand infinitesimally to the left and right of the barrier as appropriate.
Following that, we apply the diffusion equation (1) to remove the time derivative. The boundary condition in equation (3) allows us to express everything in terms of U| L , (∂U/∂ x) L , and (∂ 2 U/∂ x 2 ) L as well as p t,L→R and δ x L . Since we did not specify any constraints regarding the values of D L and D R or their relationship, p t is equivalent to p t,i→ j in equation (5) given the diffusivities D L = D i and D R = D j . The final expression is: It also suggests that the model in equation (5) as derived in 25 relies on ignoring the terms associated with U and ∂ 2 U/∂ x 2 . Since we cannot make statements about the magnitude of these terms, the former omission is only valid provided p t → 0 or D L /D R → 1, while the latter requires δ x → 0 as well. These relations are in line with observations of the model's behaviour throughout this work.

B Hybrid model derivation
First we consider the case of the membrane placed on the "slow" side, where the diffusion coefficient D 2 is lower. If the walker attempts a crossing from the slow side, it first transits through the membrane with probability P b 2 or is reflected with 1 − P b 2 . It always transits through the step change in diffusion coefficient, because it enters a region of higher diffusion coefficient. If the walker originates on the fast side, it first needs to pass the step change in diffusion coefficient with probability P d . Subsequently the membrane on the slow side either permits transit with P b 2 , which terminates the interaction, or reflects it with 1 − P b 2 , which causes the walker to attempt and always succeed in transiting through the step change in diffusion coefficient. This can be summarised with the two probabilities An alternative to the above model is to place the membrane on the "fast" side, where the diffusion coefficient D 1 is higher. This scenario results in walkers reflecting between the two interfaces continuously until eventually exiting on either side. The diagram in figure 10 illustrates this. By summing the probabilities of all different exit conditions, we obtain Here, the infinite sums only converge if |(1 − P b 1 )(1 − P d )| < 1, which is the case given that all P ≤ 1. The probabilities without power of i − 1 can be removed from the summation and, knowing that we obtain the following expressions for ...

18/22
To proof that both methods are equivalent we need to show that the following equality is valid From 25 we know that where a is a constant parameter as it only depends on κ and δt. The membrane probabilities from either side (P b 1 and P b 2 ) can then be simplified to where b and c are 1 √ D 1 and 1 √ D 2 respectively. The probability P D given by Equation (6) can be also rewritten as P D = b c .
If we consider all these simplifications, the right hand side of equation (24) is reduced to the following validating the equality in equation (24) and showing that both methods give the same probability.

C Flux analysis
Here, we further analyse the steady-state case and offer an explanation for the differences in the reference models' behaviour observed in section 3.2 and elsewhere in this work. Consider the control volume around a barrier bounded by the maximum step length δ x on either side as illustrated in figure 11. Note the difference between the running variable x (and its differential dx) and the finite step length δ x. The net flux J = −D∂U/∂ x can be split into its individual components, i.e. the amount of walkers/concentration crossing from left to right or right to left, whose magnitudes are given by: The integral bounds on the two sides cover the farthest that a walker can be located away from the boundary on either side in order for its step δ x to interact with the barrier. While the concentration density U(x) in general depends on position x, we ignore this dependency here because we assume U = constant in the domain for the steady-state case. The factor 1 /2 accounts for the fact that only half the walkers are expected to step towards the barrier. The ± signs are included to designate the directionality of the flux components J such that ∑ J = 0. For simplicity, we now omit this and use magnitudes only. We also cancel the terms U and δt and henceforth use the probability integrals F: The steady-state solution requires the net flux ∑ J = 0, i.e. F L→R = F R→L . If not, an imbalance in U(x) will develop. If the diffusivities on both sides are equal (D L = D R = D), the models considered in this work reduce to p t,L→R = p t,R→L . Furthermore, by definition δ x L = δ x R and hence ∑ J = 0. Next, we consider the different models for the non-trivial case D L = D R . Figure 11. A graphical explanation of the flux analysis around two compartments (D L > D R ) separated by a permeable barrier with asymmetric transit probabilities p t,L→R and p t,R→L . The control volume is bounded by the step lengths δ x on either side. Walkers located inside this volume at a point with concentration/density U(x) add to the flux component J with their weighted contribution ±U(x)p t (x), provided they step towards the barrier (accounted for by the factor 1 /2).

C.1 Interface model
In the case of an infinitely permeable membrane, we consider the interface model 29 . The probability of transit for interfaces with discontinuous diffusivity D (ignoring permeability κ) is given by For this model, the probabilities p t do not depend on x but only on D. As a result, we can remove p t from the integral and simplify the result: To remove the awkward min-operator for p t , consider the two distinct cases in table 1 separately. For both combinations of p t values it is clear that F L→R = F R→L . In fact, F = √ 2D min δt which we also observe in numerical simulations with κ = ∞.  (30) and their corresponding probabilities of transit p t and resulting fluxes F.

C.2 Membrane model
The authors in 25 introduce two probabilities of transit p t , going from L to R and from R to L: Each p t depends on the diffusivity D of the original compartment and the distance s from the step origin to the barrier. The latter is equal to ∓x since x b = 0, such that s L ∈ [−δ x L , 0] and s R ∈ [0, δ x R ] respectively. Note that the simplification on the right-hand sides of equation (32) only applies if 2κs D. We will analyse both cases below.

C.2.1 Larger time step
Using the full, non-approximated term for p t from equation (32), we obtain the following: The integrals can be solved to give: Since F L→R and F R→L depend on D L and D R respectively, and because F cannot be reduced further, we conclude that F L→R = F R→L and thus ∑ J = 0 as observed in section 3.

C.2.2 Sufficiently small time step
For a sufficiently small time step δt, we can solve for the simplified expression of p t : Hence, F L→R = F R→L = 2κδt (we confirm this in simulations with small δt) and the model will retain the steady-state solution.
Recall that the simplification above requires that 2κδ s D. This can be recast to 1 2κ and we now observe that the threshold at which a significant error in the fluxes occurs increases with κ and δt and decreases with D. This is in line with observations in figure 5.

C.3 Asymmetric interface reflection
As originally stated in 11 , the interface reflection condition must be respected even as κ → ∞. This is satisfied intrinsically in equation (30). Following from equation (32), the ratio of probabilities approaches unity in the limit of κ → ∞ instead of D R /D L . This might explain why the model appears to break down as the time step restriction is exceeded either by increasing δt or κ as done in figure 5. Hybrid model Reference model Figure 12. Convergence of random walk simulations of the steady-state in figure 3 using both transit models for varying number of walkers N p and time steps δt. The total simulated time is t = 100 ms and walkers are initially seeded uniformly in the domain (D L = 2.5 µm 2 /ms, D R = 0.5 µm 2 /ms, κ = 0.05 µm/ms). Figure 12 shows the histograms of walker positions for both transit models considered in this work, equations (5) and (6). We consider three time steps and increase the number of walkers until the random fluctuations are of an acceptable level (compare section 3.2).

D Convergence
.